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Avian pathogenic Escherichia coh (APEC) are responsible for heavy economic losses in poultry industry. 
Here we investigate DNA methylome of spleen and identify functional DNA methylation changes related to 
host response to APEC among groups of non-challenged chickens (NC), challenged with mild (MD) and 
severe pathology (SV). DNA methylation was enriched in the gene bodies and repeats. Promoter and CGIs 
are hypomethylated. Integration analysis revealed 22, 87, and 9 genes exhibiting inversely changed DNA 
methylation and gene expression in NC vs. MD, NC vs. SV, and MD vs. SV, respectively. IL8, IL2RB, and 
ILIRAPLI were included. Gene network analysis suggested that besides inflammatory response, other 
networks and pathways such as organismal injury and abnormalities, cell signaling and molecular transport, 
are probably related to host response to APEC infection. Moreover, methylation changes in cell cycle 
processes might contribute to the lesion phenotype differences between MD and SV. 

Avian colibacillosis is one of the most frequent infectious diseases characterized by multiple organ lesions 
typically with airsacculitis, pericarditis, perihepatitis and peritonitis''^. It is responsible for high morbidity, 
mortality and product contamination in the poultry industry, which induce significant economic losses 
and pose severe threat to human health'''\ Avian colibacillosis is caused by avian pathogenic Escherichia coli 
(APEC) infection^. APEC, a gram-negative bacterium, belongs to the extraintestinal pathogenic group of 
Escherichia coli. Although a large variety of APEC serogroups have been identified so far, the most common 
and widespread serogroups are Ol, 02 and 078^. In addition to the control of environmental factors such as 
ventilation and humidity, antibiotherapy and vaccination are the main methods for the prevention of APEC 
infection'*. However, previous research demonstrated that vaccines against heterologous APEC strains were not 
fully effective and APEC strains were frequently resistant to a wide range of antibiotics''''''''. Moreover, there is 
increasing consumer pressure to substantially decrease the use of antibiotics in food-producing animals. 
Therefore, the selection of genetically disease resistant poultry population is becoming a topic of great interest 
in the control of APEC infection. A better understanding of the host response and resistance mechanisms against 
APEC should be of great value in developing better strategies to further prevent and control APEC infection. 

Although a large number of studies in surveying the host-pathogen interactions have been conducted, most of 
them are focused on the bacterium itself to investigate virulence genes or factors involved in pathogenesis of 
APEC infection' In recent years, there were several research reported on the host transcriptomes response to 
APEC and considerable changes in gene expression were found before and after infection'^ However, so far 
little is done about the epigenetic mechanisms of host response and resistance to APEC infection. DNA methyla- 
tion is one of the main epigenetic modifications in eukaryotes. In multicellular eukaryotes, methylation occurs 
primarily at the 5-C position of cytosine within CpG dinucleotides'"". Previous studies have shown that cytosine 
DNA methylation played a central role in many biological processes like gene expression regulation, X chro- 
mosome inactivation, genomic imprinting, cancer and disease development"""^^. With the development of tech- 
nologies, a new valuable approach named Methylated DNA immunoprecipitation-sequencing (MeDIP-seq) has 
recently been extensively applied to the genome methylation studies in various species^^"^^. Moreover, MeDIP-seq 
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Table 1 | Data generated by MeDIP-Seq for each sample 

Total Unique Percentage of mapped Percentage of unique 

Sample Total number of reads Total Mapped Reads Mapped Reads reads in total reads mapped reads 

NC 36,734,694 27,087,040 17,205,739 73.74% 46.84% 

MD 36,734,694 27,062,238 17,460,699 73.67% 47.53% 

SV 36,734,694 27,375,496 17,938,237 74.52% 48.83% 



NC, MD, and SV indicated the group of non-challenged, challenged-mild pathology, and challenged-severe pathology, respectively. 



was shown to be highly efficient and reliable for methylome analysis 
with low DNA concentrations, though the resolution (about 200 bp) 
was lower than that of whole genome bisulfite sequencing^*""^'^. In 
birds, with the use of MeDIP-seq, the DNA methylation landscape 
was firstly reported in the liver and muscle tissues from the red jungle 
fowl and avian broiler™. However, their study lacked integrated ana- 
lysis with gene expression to find the potential functional genes of a 
certain trait. 

The aim of the current study was to investigate the global DNA 
methylation profiles in chicken spleen and to identify potentially 
functional DNA methylation changes related to host response and 
resistance to APEC. Here we firstly characterized the chicken 
spleen methylomes by MeDIP-seq using lUumina Hiseq 2000. 
Then, in order to validate the results of MeDIP-seq, eight regions 
of the genome were selected randomly to perform bisulfite 
sequencing experiments. Subsequently, in order to identify poten- 
tially functionally relevant DNA methylation changes, we inte- 
grated the gene expression data and DNA methylation profiles 
in each contrast of the three groups including birds with non- 
challenged (NC), challenged with mild (MD) and with severe 
pathology (SV). 

Results 

Assemble and blast analysis of MeDIP-seq reads. In this study, 
spleens of three males were used to produce one pooled DNA 
sample for each group of NC, MD and SV. A total of 36,734,694 
raw reads were generated for each of the three groups, of which more 
than 73% of the reads could be mapped and about 47% of the reads 
could be uniquely mapped to the reference chicken genome 
(Table 1). The uniquely mapping reads of NC, MD, and SV 
covered 26.6%, 25.6%, and 25.3% of the genome, respectively. 

The distribution of MeDIP-seq reads in chicken chromosomes 
(GGAl-28 and the Z chromosome) for each sample was analyzed 
and uniquely mapping reads were found in most chromosomal 
regions except for some gaps. However, in a long region of GGA17 
(from 3,180,001 to 11,182,526 bp), no unique-mapped reads but 
multi-mapped reads could be detected (Supplementary Fig. SI). 

To ascertain read distribution in different components of the gen- 
ome, we classified the uniquely mapped reads into four types: the 
reads uniquely mapped to CpG islands (CGIs), gene bodies, repeats, 
and others. Here, the criteria for CGIs identification were as follows: 
length exceeding 200 bp, GC content greater than 50%, and 
observed-to-expected CpG ratio greater than 0.6. Gene body was 
defined as the region from transcript starting site to transcript ending 
site. We found that the uniquely mapped MeDIP-seq reads were 



mainly present in the gene body regions, whereas less than 7.5% of 
the reads mapped to the CGIs (Supplementary Fig. S2). About 24% of 
the uniquely mapped reads in chicken belonged to repeat elements. 

Validation of MeDIP-seq data by bisulfite sequencing. In order to 
confirm quality of the MeDIP-seq results, three regions showing high 
methylation and one region showing low methylation were selected 
randomly to perform bisulfite sequencing in the three groups. The 
results of bisulfite sequencing were consistent with our MeDIP-seq 
data (Supplementary Fig. S3-S5). 

Global DNA methylation profUes in the chicken spleen. To assess 
overall methylation pattern in the chicken genome, we divided 
methylated regions into different components including promoters, 
gene bodies (5' UTR, 3' UTR, exon and intron), intergenic regions, 
CGI, and repeats. Here, the 2 Kb region upstream of the transcription 
start sites (TSS) were defined as the proximal promoter, while 
sequences between the 3' and 5'UTRs of the genes were defined as 
intergenic regions. In total, there were 55,978 methylated peaks in 
NC, 54,594 peaks in MD and 63,100 peaks in SV (Table 2). Peak 
distribution analysis showed that most of the peaks fell into the 
intergenic regions, 1.14% to 1.70% of the peaks were located at the 
5'UTR of genes, 2.01% to 3.13% of the peaks were at the 3' UTR, 
15.53%) to 23.61%) of the peaks were in the exons, 22.47%) to 39.44% of 
the peaks were in the intron regions, about 1 1.54% of the peaks were 
in the CGIs and more than 14.57% of the peaks were located in the 
repeats in the three groups. In addition, we investigated the relative 
methylation of each class by calculating the ratio of peaks located in 
that region to the total area of that region. Consequently, differential 
methylation levels were observed in different classes. The average 
methylation of promoters was found to be the lowest among all the 
classes and repeats exhibited a relative high level of methylation. The 
average methylation of CGIs was found to be relatively lowly 
methylated and gene bodies showed a significantly higher (P < 
0.05) level of DNA methylation than intergenic regions. Within the 
gene body, introns showed significantly higher (P < 0.05) 
methylation than UTRs and exons (Fig. lA). 

DNA methylation in promoter and gene body. On plotting 
the methylation density, low methylation levels were found in the 
proximal promoter regions and there was a sharp increase at the gene 
body regions, which were stayed at a plateau untU the transcription 
termination site (TTS) (Fig. IB). In the 2 Kb region downstream of 
TTS, we observed a clear trend for DNA methylation levels to 
increase gradually. 



Table 2 | Peak distribution in different components of the chicken genome 

Sample Total peak number promoter 5'UTR Exon Intron 3'UTR Intergenic CGIs Repeats 

NC 55978 2899 952 13219 22004 1752 38684 6598 8158 

MD 54594 1856 622 8479 12269 1097 37370 6701 8061 

SV 63100 3151 964 14659 24887 1916 43594 6668 9400 



NC, MD, and SV indicated the group of non-challenged, challenged-mild pathology, and challenged-severe pathology, respectively. 
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Figure 1 | Methylation distribution in different genomic regions. (A) . Global patterns of methylation in different genomic regions. Methylation density 
within promoter, gene body and intergenic regions was calculated by dividing the peak numbers in that region by the total area of that region. 
(B) . Methylation distribution of gene body and flanking regions. Methylation density was calculated by the ratio of methylation peak count vs number of 
data points. The gene body region referred to the region. In gene body (from TSS to TTS), each gene was split into 40 equal windows and the methylation 
density was calculated for each window. In upstream and downstream 2 kb regions, the regions were split into 20 non-overlap windows and the 
methylation density was calculated for each window. 



DNA methylation in CGIs and repeat elements. CGIs having 
methylation peaks were termed as methylated islands while the 
rest were designated unmethylated. Of the 33,915 CGIs reported in 
the chicken genome, about 19.5% (n = 6,598) were methylated in NC 
spleen, 19.8% (n = 6,701) in MD spleen, and 19.7% (n = 6,668) in SV 
spleen (Supplementary Table SI). Among these methylated CGIs, 
most were located in the intergenic regions. Within the gene body, 
exons showed higher proportion of methylated CGIs than UTRs and 
introns. And less than 2% were located in the 5'UTR region. 
Subsequently, we categorized both methylated and unmethylated 
CGIs on the basis of their sizes and calculated the CGI numbers in 
each class. About 29% of methylated CGIs were distributed in the size 
range of 200-300 bp. On the whole, the number of methylated CGIs 
decreased with the increase of size (Supplementary Fig. S6). 
Methylated CGIs were enriched in exons, while unmethylated ones 
were mainly present in introns. On the other hand, differential 
methylation was observed in different repeat types of the chicken 
genome. The predominant type of interspersed repeats was chicken 
repeat 1 - long interspersed nucleotide element (LINE/CRl), which 
accounted for more than 48% of the total methylated repeat 
sequences (Supplementary Table S2). In addition, a relatively high 
proportion (about 20%) of the methylated repeat sequences was 
observed in endogenous retrovirus like elements - long terminal 
repeat (LTR/ERVL). 

Identification of differentially methylated regions among the 
three groups. In our study, a total of 14,396 methylated genes 
were identified in the three groups, including 9,597 genes observed 
in all of the three groups. Of these, the methylation gene numbers in 
the NC, MD, and SV groups were 12,861, 10,292, and 13,418, 
respectively (Fig. 2A). Here, methylated genes were defined as 
genes overlapped with peak summits in promoter and gene body 
regions. Consequently, a total of 4,684 genes showing differential 
methylation between any two-way comparison of the groups of 
NC, MD and SV (coverage changes was more than two folds; p 
value < 0.01) were found, including 1,893 differently methylated 
genes between NC and MD (NC vs. MD), 2,970 between NC and 
SV (NC vs. SV), and 2,418 between MD and SV (MD vs. SV) 
(Fig. 2B). Moreover, 924 differentially methylated genes were 
found in both NC vs. MD and NC vs. SV, 875 in NC vs. MD and 
MD vs. SV, as well as 1,29 1 in NC vs. SV and MD vs. SV. Of these, 493 
genes were in all of NC vs. MD, NC vs. SV, and MD vs. SV. 

Identification of potentially functional relevant DNA methylation 
changes. DNA methylation in both promoters and gene bodies was 



proved to be associated with gene expression^''. Here we used the 
average normalized depth of reads as the measurement of 
methylation level for each gene. We found that gene expression 
level was negatively correlated with DNA methylation in the 
chicken genome of the three groups (Fig. 3). 

In order to identify potentially functionally relevant DNA methy- 
lation changes, we integrated the gene expression data and DNA 
methylation profiles in each contrast generated from the three groups 
of NC, MD, and SV. Differentially methylated genes with differences 
in the promoter or the gene body regions were selected to investigate 
concomitant expression changes in each contrast. Using a FDR cutoff 
of 0.001 and a filter of mean twofold change, 507 genes were found to 
be differentially expressed in NC vs. MD, 838 in NC vs. SV, and 138 
in MD vs. SV. The number of genes exhibiting coordinately changed 
DNA methylation and gene expression was described in Table 3 for 
each contrast. In the NC vs. MD contrast, a total of 11 genes were 
identified to be significandy up-methylated and down-regulated, in 
which a promoter-methylated gene and 11 body- methylated genes 
were included. Moreover, there were 11 genes significantly down- 
methylated and up-regulated, in which a promoter methylated gene 
named TXNL4B was included (Supplementary Table S3). In the NC 
vs. SV contrast, a total of 76 genes were significantly up- methylated 
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Figure 2 | Methylated genes and differentially methylated genes among 

NC, MD, and SV. (A). Methylated genes that were unique or shared among 
three groups of NC, MD, and SV. (B). Differentially methylated genes that 
were unique or shared among three groups of NC, MD, and SV. Numbers 
in each section of the figure referred to the numbers of methylated (for A) 
or differently methylated genes (for B). NC, MD, and SV indicated the 
group of non-challenged, chaUenged-mild pathology, and chaUenged- 
severe pathology, respectively. 
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Figure 3 | The correlation between DNA methylation and gene expression. X-axis is the methylation level measured by the average normalized depth of 
reads. Here methylations located in the promoter or the genebody region of genes (more than half of the sequences were overlapped with peak summits) 
were used for integration analysis. Y-axis is the gene expression level measured with RPKM value (Reads Per kb per Million reads). In general, gene 
expression level was negatively correlated with DNA methylation level. 



and down-regulated, of which 4 genes {ACAP3, ADAMTSl, CADPS2, 
CENPJ) were promoter-methylated and 73 were body-methylated. 
On the other hand, 11 genes were significandy down-methylated 
and up-regulated, of which 3 genes were promoter-methylated 
(IL8, LOC417973, TXNL4B) and 8 were body-methylated (Supple- 
mentary Table S3). Of these, five genes, PLEKHA7, ETVl, TXNL4B, 
SLC2A9, AUTS2, and ClOorfll, were found in both of the two com- 
parisons between non-chaUenged and challenged pathology groups. 
In the MD vs. SV contrast, 8 genes were statistically significantly up- 
methylated and down-regulated, while I gene was significantly down- 
methylated and up-regulated (Supplementary Table S3). All these 
genes were body-methylated. In addition to those genes showing 
inverse relationships between DNA methylation and gene expression 
changes, 7, 30, and 7 genes were found to be up-methylated but up- 
regulated, while 24, 12, and 4 genes were found to be down-methy- 
lated and down- regulated in NC vs. MD, NC vs. SV, and MD vs. SV, 
respectively (Supplementary Table S4). 



Importantly, three well-known cytokine genes, interleukin 8 {IL8), 
interleukin 2 receptor, beta {IL2RB), and interleukin 1 receptor 
accessory protein-like 1 (ILIRAPLI), exhibited significant DNA 
methylation changes and significant inverse gene expression changes 
between the NC and SV groups (Table 4). The IL8 was methylated at 
a significantly lower level (p < 0.01) and expressed at a significantly 
higher level (FDR < 0.001) in SV compared with NC. In contrast, 
both IL2RB and ILIRAPLI had a significantly higher level of DNA 
methylation (p < 0.01) and significantly lower level of expression 
(FDR < 0.001) in SV than NC. 

Subsequently, we used IPA to investigate which gene networks 
might be affected by these genes showing significant methylation 
changes in conjunction with significant inverse gene expression 
changes. The top gene network identified in the NC vs. MD contrast 
involved cellular movement, inflammatory disease and inflammat- 
ory response (Fig. 4A). Genes of ERKl/2 and ILIB were central of this 
network. A total of 1 1 genes showing inverse DNA methylation and 



Table 3 Number of genes showing both changed DNA methylation and 
Contrast Location 


gene expression in each contrast 
up- regulated 


down-regulated 


NC vs. MD 


promoter 


up-methylated 


1 


1 






down-methylated 


1 


1 




gene body 


up-methylated 


6 


11 






down-methylated 


10 


23 


NC vs. SV 


promoter 


up-methylated 


8 


4 






down-methylated 


3 


2 




gene body 


up-methylated 


22 


73 






down-methylated 


8 


10 


MD vs. SV 


promoter 


up-methylated 


0 


0 






down-methylated 


0 


0 




gene body 


up-methylated 


7 


8 






down-methylated 


1 


4 


NC vs. MD indicated the comparison between non-challenged and challenged-mild pathology groups; NC vs. SV indicated the comparison between non-challenged and challenged-severe pathology 
groups; MD vs. SV indicated the comparison between challenged-mild pathology and challenged-severe pathology groups, tocation indicated the DNA methylation differences occurred in the promoter or 
the gene body regions. Up-methylated meant that within the same peak region, there were greater reads in the second group than the first group, while down-methylated meant there were greater reads in 
the hrst group than the second group (p value < 0.01 ). Up-regulated indicated there was higher gene expression level in the second group than the first group, whereas down-regulated meant there was 
higher gene expression level in the hrst group than the second group [FDR < 0.001 ). 
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Table 4 


Comparisons of gene expression and DNA methylation between NC and SV for IL8, IL2RB, an 


d ILIRAPLI 




Gene 


Expression 




Methylation 




IL8 

IL2RB 

ILIRAPLI 


RPKM (NC) RPKM (SV) Change 
6.71 14.07 1.06 
20.52 9.67 -1.09 
1.60 0.14 -3.53 


FDR Location 
1.71E-05 Promoter 
2.98E-07 Gene body 
5.38E-05 Gene body 


Reads (NC) 
33 
26 
33 
23 


Reads (SV) 
9 
55 
74 
54 


P-value 
1 .64E-04 
2.87E-03 
2.73E-04 
1 .03E-03 


NC and SV indicated the group of non<hallenged and chollenged-severe pathology, respectively. Here change was calculated by log2 ratio (SV/NC). Location referred that the 
appeared in the promoter or gene body regions of the gene. Reads referred to normalized reads in each incorporated interval of peaks. 


differentially methylation 




Figure 4 | Identification of gene networks in the contrasts of NC vs. MD, NC vs. SV, and MD vs. SV. (A) : A gene network identified in NC vs. MD. 
(B) and (C): The top two gene networks in NC vs. SV. (D): A gene network identified in MD vs. SV. Gene network was identified through integrative 
analysis of significant DNA methylation changes in conjunction with significant inverse gene expression changes in each of the three groups. Genes 
exhibiting down-methylated and up-regulated were shown in red color, while those exhibiting up-methylated and down-regulated were shown in green. 
The intensity of red and green colors indicated the degree of up- and down-regulated, respectively. Solid lines indicated direct interaction, whereas dashed 
lines indicated indirect interaction. 
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gene expression changes were involved in this network, including 
AOAH, ETVl, ITGA9, LARP6, PLCXDl, PLEKHA7, PTPRZl, 
SLC2A9, SLC8A1, SPIl, and TSPAN12. In the NC vs. SV contrast, 
the top two gene networks identified involved cancer, organismal 
injury and abnormalities, on the one hand, and cell signaling, 
molecular transport, on the other (Fig. 4B, 4C). Of those genes show- 
ing inverse relationships between DNA methylation and gene 
expression changes, a total of 18 genes {ADAMTSl, ATP8A2, 
AUTS2, CDS A, CENP}, DOCKS, EDAR, FREMl, FSHR, KIAA0226, 
KIT, MICAL2, NEK3, PDEIC, SLC2A9, ST3GAL1, TACRl, and 
TPMl) were involved in the first network, while 16 genes {BPI, 
CACNAIH, CNGA3, ETVl, F3, GRIN2B, IL8, ITPRl, KDMS, 
LRP5, NFAMl, NFATCl, PRKDl, SLC8A1, TMEFF2, and WASF3) 
were involved in the second network. Importantly, the cytokine IL8 
was one of the central molecules in the second network. On the other 
hand, one gene network involved in cell cycle was identified in the 
MD vs. SV contrast and all genes we obtained (CCT5, CENP], 
CAMK2G, DNAAF2, LAP3, MYSMl, PTCHD2, SCUBEl, and 
VIPRl) were involved in this network (Fig. 4D). 

Discussion 

MeDIP-seq peaks are the methylation enriched regions. Peak scan- 
ning and analysis were important to survey the methylome pattern. 
In the present study, we assessed the chicken spleen methylation 
pattern through the analysis of the genomic distribution of peaks 
and found that DNA methylation was enriched in the gene body 
regions and repeat elements. These findings confirmed the results 
of previous study conducted by Li et al.^° using the chicken liver and 
muscle tissues by the MeDIP-Seq method. In mammals, both pro- 
moter and CGIs were observed to be hypomethylated''^ '^, whereas 
the methylation in gene bodies was found to be higher than that in 
intergenic regions'"^-'". AU of these results were consistent with the 
results of the current study. However, in contrast to previous studies 
in animals^^'"*^'^', we did not observe that exons exhibited higher 
methylation level than introns in chickens. 

The repressive effects of promoter-associated methylation on gene 
expression have been well demonstrated'''. Earlier reports revealed 
that the methylation pattern of promoter regions showed a V shaped 
curve indicative of low methylation levels at the TSS in human, mice, 
as well as rat, whereas genes in the puffer fish exhibited a prominent 
dip in methylation just upstream of the TSS^^''^'"" ''. In many plants 
such as Arabidopsis and rice, CG methylation exhibited a character- 
istic peak in the body of protein-coding genes''. In the present study, 
DNA methylation was lowest at about 400 bp upstream of the TSS, 
then increased sharply and stayed at a plateau in the gene body 
regions, which were in accordance with the earlier investigation for 
chicken'". Previous studies demonstrated that DNA methylation in 
the gene body regions might alter transcription elongation efficiency 
and regulate cell context-specific alternative promoters in gene 
bodies"""*. The relatively high methylation level of the gene body 
in the chicken suggested that the methylation of the gene body 
regions played an important role in regulating gene expression. 

In contrast to the 40-50% density observed in mammalian gen- 
omes, the density of interspersed repeats in chicken genome is less 
than 9%" The predominant interspersed repeat element is the 
CRl long interspersed nucleotide element (LINE) and it accounts 
for over 80% of all interspersed repeats in the chicken genome'^. 
Based on our methylome analysis in chicken spleens, we identified 
LINE/CRl to be the predominant target of DNA methylation, 
accounting for more than 48% of the total methylated repeats and 
these data were consistent with the recent study of the chicken DNA 
methylome'". High-throughput DNA sequencing revealed similar 
numbers of CGIs in humans and mice: 25,495 and 23,02 1 per haploid 
genome, respectively In this study, a total number of 33,915 CGIs 
were identified in the chicken genome. Previous reports demon- 
strated that most of the CGIs were unmethylated and our results 



for chicken confirmed these earlier conclusions'" ''''. The methylation 
of the 5'UTR region usually led to the suppression of gene express- 
ion. Our data showed that less than 2% of methylated CGIs were 
located in this region in the chicken genome. Although slightly lower, 
these results were consistent with the observations reported by pre- 
vious study'". It has been shown that the methylation of promoter 
CGIs would induce transcriptional silencing of their downstream 
genes by changing chromatin structures and blocking transcription 
initiation" ''^. Recently, many studies have uncovered a large class of 
CGIs named orphan CGIs, which were remote from annotated pro- 
moters and could be categorized as intragenic or intergenic CGIs, 
and confirmed that these CGIs had the characteristics of functional 
promoters'"'''^'"'. Moreover, further analysis suggested that the 
majority of methylated CGIs were located in intragenic and inter- 
genic regions'"'''^. In the present study, we also found that methylated 
CGIs in chicken were mainly present in the intergenic regions fol- 
lowing by exons. On the other hand, a recent study conducted in rat 
showed that a large proportion of methylated CGI were located in the 
size range of 200-300 bp and the number of CGIs decreased with 
increase in the size of the islands, which were in accordance with the 
results of the current study. 

The whole chicken genome sequenced in 2004 is 1,063 Mb in total 
length and contains 20,000-23,000 predicted genes". In our study, a 
total of 14,396 methylated genes were identified and about two- 
thirds of them were found in all of the three groups. These findings 
suggested that DNA methylation is a common feature of the chicken 
genome. Moreover, our results demonstrated that DNA methylation 
level was negatively correlated with gene expression level. The spleen 
is an important immunological organ in chicken and it has been 
successfully utilized for studies about host immune response to dis- 
ease"-''*'"'. Therefore, here we used spleen tissues to detect DNA 
methylation differences after challenged with APEC. Based on the 
comparison analysis among the NC, MD and SV groups, we iden- 
tified 4,684 genes showing significant differential DNA methylation 
in spleens. These genes were useful for the following detection of 
potential genes affecting host response and resistance to APEC. 
Sandford et al." studied the spleen transcriptome response to 
APEC infection in chicken and found very little expression difference 
between mildly infected and non-infected groups on either day 1 or 
day 5 post infection. Nie et al.''' also found fewer significantly differ- 
entially expressed unigenes in the chicken spleen transcriptome of 
NC vs. MD than that of NC vs. SV. In the current study, methylome 
comparison analysis among different treatment groups indicated 
that the birds with severe pathology would produce larger DNA 
methylation changes. 

By integrating DNA methylation and mRNA expression data, a 
number of genes exhibiting coordinately changed DNA methylation 
and gene expression were identified in this study. Of these, 22, 87, 
and 9 genes showed significantly inversely correlation between 
changes of DNA methylation and changes of expression in the NC 
vs. MD, NC vs. SV, and MD vs. SV contrasts, respectively. These 
DNA methylation changes might be functional relevant changes and 
those genes might be useful for future epigenetic studies about host 
response and resistant to APEC infection. Some of those genes, such 
as FABP3 and IL8, were investigated in previous studies of chicken 
response to virus infection'''^™. However, for most of them, so far little 
was known about DNA methylation-based deregulation in chickens 
challenged by pathogen. 

Our study revealed that genes with inversely changed DNA 
methylation and gene expression in different contrasts enriched in 
different biological pathways. These results indicated that the pos- 
sible mechanism of host response to pathogen challenge might be 
distinct in a particular pathology state. The network related to 
inflammatory and cellular movement appeared in the contrast com- 
pared between birds with mild infection and control, and 1 1 genes 
were implicated in this network. Inflammation is the first response to 
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infection and it is an essential component of the organism's defense 
mechanism against pathogens and tissue damage^'. The purpose of 
inflammatory response is to control or eradicate the infection and 
heal the damage. The network we found in NC vs. MD suggested that 
DNA methylation changes might play important roles in modulating 
the host immune/inflammatory response to APEC via the regulation 
of gene expression. On the other hand, in the contrast compared 
between birds with severe lesions and control, top gene network 
associated with organismal injury and abnormalities was observed 
and 18 genes were implicated in it, suggesting that protective res- 
ponses to organismal injury were activated after severe infection and 
methylation changes in those genes might have functional conse- 
quences in organismal injury/survival. In addition, network invol- 
ving cell signaling and molecular transport in the NC vs. SV contrast 
highlighted the importance of proper signaling cascades to fight 
severe infection. In the current study, a network involving cell cycle 
was found in the comparison between the MD and SV chickens, 
which showed different response to APEC infection, and 9 of the 
genes we identified were contained in this network. These findings 
indicated that the DNA methylation changes in the host cell cycle 
processes might contribute to the lesion phenotype differences. 

Among the genes showing significant methylation changes in 
conjunction with significant inverse gene expression changes, three 
pro-inflammatory cytokine genes, IL8, IL2RB, and ILIRAPLI, were 
of particular interest. Cytokines were small proteins released by 
immune system cells and other cells, and they were an important 
part of the intercellular communication system responsible for 
immune response^^. Generally, cytokines fulfilled the biological pro- 
cesses by binding their specific receptors. Chemokines were a special 
type of cytokines that induced the migration of cells to sites of infec- 
tion or injury""^. The interleukin molecules (ILs) were an important 
class of cytokines. Some ILs, such as ILS, also belonged to chemo- 
kines^^. Recently, ILS was reported to play an important role in 
immune response to Salmonella in chicken^'. In addition, it was 
demonstrated to be one of important genes involved in chicken 
MD-resistance or -susceptibility^". In the present study, ILS was a 
central gene of the network involving cell signaling and molecular 
transport in NC vs. SV. Birds with severe pathology had significantly 
lower levels of DNA methylation and significantly higher levels of 
expression of ILS than birds with non-challenged. Moreover, the 
altered DNA methylation was occurred in the promoter region of 
the ILS gene. These indicated that the DNA methylation change of 
ILS might play a critical role in modulating the host response and 
resistance to APEC. ILIRAPLI, similar to the interleukin 1 accessory 
proteins, was a member of the interleukin 1 receptor family. Earlier 
study has demonstrated an involvement of ILIRAPLI in the immune 
proinflammatory response to pathogen infection in pigs". On the 
other hand, IL2RB was reported to play a key role in T-ceU mediated 
immune response and its genetic polymorphisms were associated 
with various diseases^^'^''. In our study, both ILIRAPLI and IL2RB 
were significantly up-methylated and down-regulated in SV com- 
pared to NC, which supports the critical function of DNA methyla- 
tion changes of these two genes in the immune response to APEC in 
chicken. 

In conclusion, we have generated the splenic DNA methylome for 
the chicken and the methylated regions obtained from MeDIP-seq 
could be validated by bis-seq. The genome-wide DNA methylation 
profiles were compared among the NC, MD, and SV groups. By 
integrating DNA methylation and mRNA expression data, a number 
of potentially functional relevant DNA methylation changes were 
identified, of which some played important roles in the regulation 
of the expression of genes involved in the host inflammatory 
response and organismal injury/survival after APEC infection. 
Methylation changes of genes involving cell signaling and molecular 
transport were also important in fighting severe infection. Moreover, 
DNA methylation changes of genes in the cell cycle processes might 



contribute to the lesion phenotype differences between MD and SV 
chickens. 

Methods 

Ethics statement. All experiments of this study were approved by the Iowa State 
University Institutional Animal Care and Use Committee (# 11-07-6460-G). And the 
experiment was performed according to regulations and guidelines established by this 
committee. AU efforts were made to minimize suffering. 

Animals and sample collection. A total of 360 non-vaccinated commercial male 
broilers were reared in cages with ad libitum access to water and food. All chickens 
were subj ectedtoa22:2hour light ; dark cycle for the first 1 5 days and then changed 
to a 16 : 8 hour cycle. At 4 weeks of age, birds were divided into two groups. One group 
(n — 240) was challenged with 0.1 ml {containing 108 colony forming units) APEC 
Ol. The others (n — 120) were assigned to the control group, in which birds were 
non-challenged (NC) but treated with 0.1 ml phosphate buffered saline (PBS). At one 
day post challenge, birds were sacrificed and necropsied to determine their lesion 
scores as described by previous study^"*. Birds with scores ranging from 0 to 2 were 
designated as mild pathology, while those with scores ranging from 4 to 7 were 
regarded as severe pathology. The average lesion scores for NC, JVID, and SV were 0.00 
± 0.00, 0.50 ± 0.58, and 5.25 ± 1.26, respectively. In each group, the whole spleens 
from 3 individuals were collected and placed in RNAlater immediately, and then 
stored at — 80 C untU DNA extraction. 

DNA preparation and MeDIP-seq. Genomic DNA was extracted using Qiagen 
DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA, USA) following the 
manufacturer's recommendation and then treated with 10 fig/ml RNase enzyme 
(Qiagen, Valencia, CA, USA) to degrade RNA. The quality and purity of DNA were 
determined by agarose gel electrophoresis and spectrophotometer. DNA quantity 
was measured with Quant-iT dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). 
DNA from three individuals within each group was pooled in equal amounts to 
produce a mixed sample per group. Subsequently, 5 |ig mixed DNA sample was 
sonicated to generate random fragments ranging from 100 to 500 bp. Then end- 
repair, phosphorylating and A-taUing were performed with Paired-End DNA Sample 
Prep kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocol. 
After ligation with Illumina adaptors, the fragments were used for subsequent MeDIP 
enrichment with Magnetic Methylated DNA Immunoprecipitation kit (Diagenod, 
Liege, Belgium) using 1 [ig 5-methylcytosine antibody. DNA from qualifying MeDIP 
experiment was amplified to produce libraries with insert sizes between 200 and 
300 bp. Products were quantified using an Agilent 2100 Analyzer (Agilent 
Technologies, Santa Clara, CA, USA) and the Quant-iT dsDNA HS Assay Kit 
(Invitrogen, Carlsbad, CA, USA). Sequencing was carried out on the Illumina Hiseq 
2000 (Illumina, San Diego, CA, USA) following the standard protocol to generate 
paired-end 50-bp reads in a commercial company (BGI, Shenzhen, Guangdong, 
China). 

Bisulfite sequencing. One microgram of pooled DNA from each group was bisulfite- 
treated using the EpiTect Bisulfite kit (Qiagen, Valencia, CA, USA) according to the 
manufacturer's instructions. Four pairs of primers were designed with Methyl Primer 
Express Software vl.O. PI (F: TTGTTTTTGTTTTGGATGGTTA, R: 
ATTAACCACCCCCACCTAC), P2 (F: AGGGTTGTAGATTTGTATTTGGA, R: 
ACTTAATACCTTTCTCCCCATCT) and P3 (F: GTTGATTGTAGTT- 
TTTTGAAGG, R: CACATCCTCAATAATAATATCCC) were used for the 
validation of regions with high methylation, whereas P4 (F; ATGTGGTATTGAGG- 
GATATAGGT, R: TAAATAATTTCCCCCCCTATC) was used for the vahdation of 
regions with relatively low methylation (Supplementary Table S5). Semi-nested PGR 
was performed in 50-[il reaction volumes containing 50 ng template DNA, 1 |iM of 
each primer and 2.5 U of LA Taq HS (TaKaRa, Osaka, Japan). The first round and the 
second round of amplification conditions were the same as follows; 94 'C for 3 min; 
35 cycles of 94"C for 30 s, 56 C for 30 s and 72°C for 30 s; and 72^C for 5 min. PGR 
products were gel-purified with a Gel Extraction Kit (Tiangen, Beijing, China) and 
then cloned into the pMD18-T vector (Takara, Osaka, Japan). Ten clones were 
sequenced and then analyzed using ClustalW. 

RNA-Seq. The process of RNA-Seq was described in our previous study^'^. Briefly, 
total RNA of spleens used for the DNA methylation analysis was extracted with 
TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Then, three spleens in each group 
were pooled in equal amounts to generate one mixed sample and these three mixed 
samples were subsequently used for cDNA library construction with the TruseqTM 
RNA sample prep Kit (Illumina, San Diego, CA, USA) according to the 
manufacturer's instructions. After the steps of size selection, PGR amplification, 
product validation and quantification, libraries were subjected to Illumina deep 
sequencing on the Genome Analyzer IIx (Illumina, San Diego, CA, USA) in Shanghai 
Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China). 

MeDIP-Seq data analysis. Reads containing adapters, unknown or low quality bases 
were filtered out from raw data generated by Illumina sequencing. Then the retained 
reads were mapped to the chicken reference genome downloaded from Ensembl 
genome browser (ftp://ftp.ensembl.org/pub/release-63/fasta/gallus_gallus/dna/) 
using the software SOAPaligner v 2.21 (http;//soap. genomics. org.cn/)^^. Mismatches 
no more than 2 bp were allowed in the alignment and uniquely mapped reads were 
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retained for further analysis. Gene information was also obtained from the public FTP 
site of Ensembl (ftp://ftp.ensembl.org/pub/release-63/gtf/gallus_gallus/). Repeat 
annotations were downloaded from the UCSC database (http://hgdownload.cse.ucsc. 
edu/goldenPath/rn4/bigZips/chromOut.tar.gz) and reads distribution on repetitive 
elements was analyzed by the software RepeatMasker (http://www.repeatmasker.org/ 
). The CGIs were scanned by CpGPlot (https://gcg.gwdg.de/emboss/cpgplot.html) in 
this study. The peak detection was carried out with the software MACS V 1.4.2 (http:// 
liulab.dfci.harvard.edu/MACS/) Statistical analyses of methylation level differences 
in different classes was performed with least square method using SAS 8.0 software 
(SAS Institute Inc., Cary, NC, USA). All genes with peaks were employed for the 
subsequent identification of differentially methylated genes that exhibited more than 
2-fold changes in the three samples while compared to each other and p value < 0.0 1 . 
Then GO enrichment analysis was conducted for those genes exhibiting altered DNA 
methylation using the DAVID Functional Annotation Tool with a 0.05 cutoff for 
Benjamini adjusted p-value^^. 

Integration of gene expression analysis. For the raw reads obtained from RNA-Seq, 
the sequencing adaptors, reads with unknown nucleotides larger than 5%, and low 
quality bases (more than half of the bases' qualities were less than 10) were trimmed. 
Then all the clean data was mapped to the chicken genome with no more than 5 bp 
mismatches using SOAPaligner v 2.21 (http://soap.genomics.org.cn/)^^. Sequences 
uniquely mapped to the chicken genome were used for subsequent analysis. Gene 
expression level was calculated with RPKM method (Reads Per kb per Million reads)^" 
and genes with more than twofold changes and false discovery rate (FDR) ^ 0.001 
were regarded as significant differentially expressed. Differentially methylated genes 
located in the promoter or the genebody region genes (more than half of the 
sequences were overlapped with peak summits) were used to identify potentially 
functional genes affecting host response and resistance to APEC by integrating with 
mRNA expression data. Additionally, genes of which expression levels changed in 
accordance with DNA methylation changes were used for gene network and 
biological processes enrichment analysis with Ingenuity Pathways Analysis (IPA) 
(Ingenuity Systems; http://www.ingenuity.com). 
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